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Ising model with quenched random magnetic fields is examined for single Gaussian, bimodal 
and double Gaussian random field distributions by introducing an effective field approximation 
that takes into account the correlations between different spins that emerge when expanding the 
identities. Random field distribution shape dependence of the phase diagrams, magnetization and 
internal energy is investigated for a honeycomb lattice with a coordination number q — 3. The 
conditions for the occurrence of reentrant behavior and tricritical points on the system are also 
discussed in detail. 
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I. INTRODUCTION 

Ising model (T), [1] which was originally introduced as a 
model describing the phase transition properties of ferro- 
magnetic materials has been widely examined in statis- 
tical mechanics and condensed matter physics. In the 
course of time, basic prescience of this simple model 
has been improved by introducing new concepts such as 
disorder effects on the critical behavior of the systems 
in question. Ising model in a quenched random field 
(RFIM) which has been studied over three decades is an 
example of this situation. The model which is actually 
based on the local fields acting on the lattice sites which 
are taken to be random according to a given probabil- 
ity distribution was introduced for the first time by Imry 
and Ma Q. Subsequently, a great many of theoretical 
and experimental works have paid attention on the RFIM 
and quite noteworthy results have been obtained. For in- 
stance, it has been shown that diluted antiferromagnets 
such as Fe x Zni- x F 2 $, 0, Rb 2 Co x Mgi^ x F^ 0] and 
Co x Zni- x F2^ in a uniform magnetic field just corre- 
spond to a ferromagnet in a random uniaxial magnetic 
field [1[9|. 

Following studies have been devoted to investigate the 
phase diagrams of these systems in depth and in the mean 
field level it was found that different random field dis- 
tributions lead to different phase diagrams. For exam- 
ple, using a Gaussian probability distribution Schneider 
and Pytteplj have shown that phase diagrams of the 
model exhibit only second order phase transition prop- 
erties. On the other hand, Aharony [ll[ and Matthis 
[i~2| have introduced bimodal and trimodal distributions, 
respectively and they have reported the observation of 
tricritical behavior. In order to clarify this controversial 
situation, the problem have been investigated by a va- 
riety of theoretical works such as effective field theory 
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(EFT) 113423, Monte Carlo (MCj simulations |24>|29j, 
mean field (MF) approximation [3(M35| , pair approxima- 
tion (PA) [36|, |37|, Bethe-Peierls approximation (BPA) 
HU and series expandion (SE) [3^ method. Moreover, 
recently phase transition properties of RFIM with sym- 
metric double [!o[ and triple [4lJ Gaussian random fields 
have also been studied by means of a replica method and 
a rich variety of phase diagrams have been presented. 

EFT has been widely used in the literature and some 
interesting results have been reported. For example, de 
Albuquerque et al. [H[ have studied the behavior of the 
spin-1/2 RFIM on a honeycomb lattice with a bimodal 
field distribution and they have observed the absence 
of a tricritical point. Liang and co-workers fl9| exam- 
ined spin-3/2 system with a random field on honeycomb, 
square, and simple-cubic lattices, respectively and re- 
ported the existance of tricritical point and the reentrant 
phenomena for a system with any coordination number- 
ed when the ap plie d random field is bimodal. Sarmento 
and Kaneyoshi [21] have studied the phase diagrams of a 
transverse spin-1/2 Ising model with a random field and 
found a reentrant behavior of second order. They have 
not observed any tricritical point for the system with a 
coordination number q = 3 when the applied field is bi- 
modal. Sebastianes and Figueiredo [22| have studied a 
semi-infinite simple cubic lattice with a trimodal random 
field distribution and they have observed first order tran- 
sitions on the system. 

Conventional EFT approximations include spin-spin 
correlations resulting from the usage of the Van der Waer- 
den identities and provide results that are superior to 
those obtained within the traditional MFT. However, 
these conventional EFT approximations are not sufficient 
enough to improve the results, due to the usage of a de- 
coupling approximation (DA) that neglects the correla- 
tions between different spins that emerge when expand- 
ing the identities. Therefore, taking these correlations 
into consideration will improve the results of conven- 
tional EFT approximations. In order to overcome this 
point, recently we proposed an approximation that takes 
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into account the correlations between different spins in 
the cluster of considered lattice [42Tj46| . Namely, an ad- 
vantage of the approximation method proposed by these 
studies is that no uncontrolled decoupling procedure is 
used for the higher order correlation functions. On the 
other hand, as far as we know EFT studies in the lit- 
erature dealing with RFIM are based only on discrete 
probability distributions (bimodal or trimodal). Hence, 
in this work we intent to study the phase diagrams of the 
RFIM with single Gaussian, bimodal and double Gaus- 
sian random field distributions. 

Organization of the paper is as follows: In section UU 
we briefly present the formulations. The results and dis- 
cussions are presented in section Hill and finally section 
IIVI contains our conclusions. 



II. FORMULATION 

We consider a two-dimensional lattice which has N 
identical spins arranged. We define a cluster on the lat- 
tice which consists of a central spin labeled So, and q 
perimeter spins being the nearest neighbors of the central 
spin. The cluster consists of (q+ 1) spins being indepen- 
dent from the value of S. The nearest neighbor spins are 
in an effective field produced by the outer spins, which 
can be determined by the condition that the thermal av- 
erage of the central spin is equal to that of its nearest 
neighbor spins. The Hamiltonian describing our model 
is 



H = - J E s'S-i-'EhiS; 



(i) 



<i,3> 



where the first term is a summation over the nearest 
neighbor spins with Sf = ±1 and the second term rep- 
resents the Zeeman interactions on the lattice. Random 
magnetic fields are distributed according to a given Gaus- 
sian distribution function. Present study deals with three 
kinds of Gaussian field distribution. Namely, a normal 
distribution which is defined as 
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with zero mean and width <r, a bimodal discrete distri- 
bution 

P(hi) = \ [S(hi - h ) + <5(^ + h )} , (3) 

where half of the lattice sites subject to a magnetic field 
ho and the remaining lattice sites have a field — ho, and 
a double peaked Gaussian distribution 
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In a double peaked distribution defined in equation ((4]), 
random fields ±h are distributed with equal probability 
and the form of the distribution depends on the ho and 
a parameters, where a is the width of the distribution. 

According to Callen identity [47[ for the spin- 1/2 Ising 
ferromagnetic system with the coordination number q, 
the thermal average of the spin variables at the site i is 
given by 



({fi}Sf)= {7*}tanh 



Sj 



, (5) 



where V = d/dx is a differential operator, j expresses the 
nearest neighbor sites of the central spin and {fi} can be 
any function of the Ising variables as long as it is not a 
function of the site. From equation ([5]) with = 1, the 
thermal and random-configurational averages of a cen- 
tral spin can be represented in the form for honeycomb 
(q = 3) lattice by introducing the differential operator 
technique [H, [|| 

« S o»r = //njcosh(JV)+5| S inh(JV)]\\ 

xF(x)| x=0 ,(6) 

where the inner (...) and the outer (...) r brackets repre- 
sent the thermal and configuratio- 

nal averages, respectively. The function F(x) in equation 
^ is defined by 



F{x) = dhiP(hi)taxih.[P(x + hi)}, 



(7) 

and it has been calculated by numerical integration and 
by using the distribution functions defined in equations 
((21), ([3]) and Q. By expanding the right hand side of 
equation ([6]) we get the longitudinal spin correlation as 

((S5))r - fco + 3fci((^i)) r + 3fc 2 ((,S 1> 5' 2 }} J . 

+k 3 ((S 1 S 2 S 3 )) r . (8) 

The coefficients in equation ([5]) are defined as follows 

fc = cosh(JV)F(x)| 3;= o, 

fci = cosh 2 (JV)sinh(JV)F(x)| a;= o, 

k 2 = cosh(JV)sinh 2 (JV) J F(x)| a;= o, 

k 3 = sinh 3 (JV)F(a;)U = o. (9) 

Next, the average value of the perimeter spin in the sys- 
tem can be written as follows and it is found as 

77ii = ((5 , r))r = «cosh( JV) + smh(JW))) r F(x + 7 ), 
= ai + a 2 ((S )) r . (10) 

For the sake of simplicity, the superscript z is omitted 
from the right hand sides of equations ((8|) and (fTOf . The 
coefficients in equation (fT0|) are defined as 



(11) 



ai = cosh(JV) J F(x + 7)| x =o, 
a 2 = sinh(JV)F(a; + 7)^=0. 



3 



In equation 7 = (q — 1)A is the effective field 

produced by the (q — 1) spins outside of the system 
and A is an unknown parameter to be determined self- 
consistently. Equations (JSJ and (fTU|) are the fundamental 
correlation functions of the system. When the right-hand 
side of equation (|6]) is expanded, the multispin correla- 
tion functions appear. The simplest approximation, and 
one of the most frequently adopted is to decouple these 
correlations according to 

((S?(S*f...S? » r = «s?» r «(s*) 2 )) r ... {(Sf » r , 

_ (12) 
for i ^ j ' ^ ... ^ I [50|. The main difference of the method 
used in this study from the other approximations in the 
literature emerges in comparison with any decoupling ap- 
proximation (DA) when expanding the right-hand side of 
equation ©. In other words, one advantage of the ap- 
proximation method used in this study is that no uncon- 
trolled decoupling procedure is used for the higher order 
correlation functions. 

For spin- 1/2 Ising system in a random field with q = 3, 
taking equations ([5]) and (ITU1) as a basis we derive a set of 
linear equations of the spin correlation functions which 
interact in the system. At this point, we assume that (i) 
the correlations depend only on the distance between the 
spins and (ii) the average values of a central spin and its 
nearest-neighbor spin (it is labeled as the perimeter spin) 
are equal to each other with the fact that, in the matrix 
representations of spin operator S, the spin-1/2 system 
has the property (5|) 2 = 1 . Thus, the number of linear 
equations obtained for the system with q — 3 reduces to 
six and the complete set is as follows 

((5 »r = fco + 3fci((5 1 )) r + 3fc 2 ((5 1 5 2 )) r 

+ k 3 ((S 1 S 2 S 3 )) r 

((S Si)} r = 3k 1 + (k + 3k 2 )((S 1 )) r + k 3 ((S 1 S 2 )) r 

((SoS 1 S 2 )) r = (3k 1 + k 3 )((S 1 )) r + (k + Sk 2 ){(S 1 S 2 )) r 

((S 1 )) r = a 1+ a 2 ((S Q )) r 

((SiS 2 )) r = a 1 ((S 1 )) r + a 2 ((S S 1 )) r 

((SiS 2 S 3 )) r = ai((5 1 5 2 )) r + a 2 ((5 5 1 5 2 )) r . (13) 

If equation (|13|) is written in the form of a 6 x 6 
matrix and solved in terms of the variables Xi[(i — 
1,2, ...,6)(e.g.,xi = {(S§)) r ,x 2 = ((S Si)) r , x 6 = 
{{SiS 2 S 3 )) r )} of the linear equations, all of the spin cor- 
relation functions can be easily determined as functions 
of the temperature and random field parameters a and 
ho. Since the thermal and configurational averages of the 
central spin is equal to that of its nearest-neighbor spins 
within the present method then the unknown parameter 
A can be numerically determined by the relation 

({S ))r = ((Si))r or X t = Xt. (14) 

By solving equation (I14p numerically at a given fixed set 
of Hamiltonian/random field parameters we obtain the 
parameter A. Then we use the numerical values of A to 
obtain the spin correlation functions which can be found 



from equation (fl3l) . Note that A = is always the root 
of equation (fl"4"|) corresponding to the disordered state of 
the system. The nonzero root of A in equation (fl"4"|) cor- 
responds to the long range ordered state of the system. 
Once the spin correlation functions have been evaluated 
then we can give the numerical results for the thermal 
and magnetic properties of the system. Since the effec- 
tive field 7 is very small in the vicinity of fcsTc/ J, we can 
obtain the critical temperature for the fixed set of Hamil- 
tonian parameters by solving equation (j!4[) in the limit of 
7 — > then we can construct the whole phase diagrams 
of the system. Depending on the Hamiltonian parame- 
ters, there may be two solutions (i.e. two critical temper- 
ature values satisfy the equation (O) corresponding to 
the first/second and second order phase transition points, 
respectively. We determine the type of the transition by 
looking at the temperature dependence of magnetization 
for selected values of Hamiltonian parameters. 



III. RESULTS AND DISCUSSION 

In this section, we discuss how the type of random 
field distribution effects the phase diagrams of the sys- 
tem. Also, in order to clarify the type of the transitions 
in the system, we give the temperature dependence of 
the order parameter. 



A. Phase diagrams of single Gaussian distribution 

The form of single Gaussian distribution which is de- 
fined in equation @ is governed by only one parameter 
a j J which is the width of the distribution. On the left 
panel in Fig. la, we show the phase diagram of spin-1/2 
system on a honeycomb lattice in (ksT c / J ~ a/ J) plane. 
We can clearly see that as a j J increases then the width of 
the distribution function gets wider and randomness ef- 
fect of magnetic field distribution on the system becomes 
significantly important (see the right panel in Fig. la). 
Therefore, increasing a/J value causes a decline in the 
critical temperature ksT c / J of the system. ksT c / J value 
in the absence of any randomness i.e. when a/J = is 
obtained as ksTc/J — 1.5044 and this value just corre- 
sponds to the critical temperature of the pure system. 
This value is in a good agreement with the exact result 
1.519 [5l|| . We also note that the critical temperature 
of the system reaches to zero at a/J — 1.435. Besides, 
we have not observed any reentrant/tricritical behavior 
for single Gaussian distribution. The system undergoes 
only second order phase transition as depicted in Figs, 
lb and lc. In Figs, lb and lc, we show the order param- 
eter (magnetization) and the internal energy surfaces in 
(m-k B T/J-a/J) and {U / J N - k B T / J - a / J) planes, 
respectively. As we can clearly see from Fig. lb, as the 
temperature increases then the magnetization of the sys- 
tem tends to decrease from its saturation value and falls 
to zero at a second order phase transition temperature. 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 
o/J 



(a) 




(b) 




(c) 



FIG. 1: (a) Phase diagram (left panel) of system corre- 
sponding to single Gaussian distribution in (fcsT c /J — o/J) 
plane and variation of single Gaussian distribution with ho /J 
for selected values of o/J (right panel), (b) Magnetization 
and (c) internal energy surfaces on (m — fcsT/J — o/J) and 
(U/JN — ksT / J — o/J) planes, respectively. 



Furthermore, in the presence of any random field corre- 
sponding to equation @, increasing o/J values give rise 
to a decrement in the critical temperature of the system 
and hence the ferromagnetic region gets narrower. The 
internal energy as a function of the temperature and dis- 
tribution width a / J which is depicted in Fig. lc also 



agrees with the observations mentioned above. Thus, 
we see that our results support the previously published 
works in the literature. 



B. Phase diagrams of bimodal distribution 

Next, in order to investigate the effect of the bimodal 
random fields defined in equation © on the phase dia- 
grams of the system, we show the phase diagram on the 
(A;bT c /J — ho/ J) plane and temperature dependence of 
magnetization m on the left and right hand side panels in 
Fig. 2, respectively. The right panel in Fig. 2 shows three 
different magnetization profiles. Namely, a second order 
transition occurs between paramagnetic and ferromag- 
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FIG. 2: Phase diagram in (fcsT c — ho/ J) plane (left panel) 
corresponding to bimodal field distribution obtained from 
present work and results of Ref.[l4| and Ref.[2l| and temper- 
ature dependence of magnetization (right panel) with some 
selected values of ho/ J. Solid and dotted lines correspond 
to second and first order phase transitions and solid circle 
denotes the tricritical points. 

netic phases for ho/ J — 0.9, a first order phase transition 
(dotted line on the right panel in Fig. 2) is followed by a 
second order phase transition for ho/ J = 1.08 and a reen- 
trant behavior of second order appears for ho/ J = 1.19. 
On the left panel in Fig. 2, we also compare our results 
by the other works based on EFT in the literature. In 
Refs.[l4[ and [2l| the authors reported a reentrant behav- 
ior for spin-1/2 system on a honeycomb lattice (q = 3) 
with bimodal random fields but they did not observe a 
tricritical point. On the other hand, as seen from Fig. 2 
the system exhibits a tricritical behavior as well as reen- 
trant phenomena within the present work. According 
to our calculations, reentrant phenomena and the first 
order phase transitions can be observed in the range of 
1.0 < h / J < 1.215 and 1.0 < h /J < 1.151, respec- 
tively. Furthermore, for ho/ J = (i.e. pure system 
again) the critical temperature kpT r JJ — 1.5044 can be 
compared with 1.9103 of Ref.pJ] and 2.104 of Ref.[21| 
and the exact result 1.519 of Ref . [5l| . Hence, it is ob- 
vious that the present method is superior to the other 
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EFT methods in the literature. The reason is due to the 
fact that, in contrast to the previously published works 
mentioned above there is no uncontrolled decoupling pro- 
cedure used for the higher order correlation functions 
within the present approximation. 



C. Phase diagrams of double Gaussian distribution 

Double Gaussian distribution in equation (|4]) have not 
yet been examined within the framework of EFT in the 
literature. Therefore, it would be interesting to investi- 
gate the phase diagrams of the system with random fields 
corresponding to equation Now the shape of the ran- 
dom fields is governed by two parameters h / J and a/J. 
In Fig. 3a, right panel shows three different magneti- 
zation profiles for a/J — 0.1. Namely, a second order 
transition occurs between paramagnetic and ferromag- 
netic phases for ho/ J = 1.0, a first order phase transition 
(dotted line on the right panel in Fig. 3a) is followed by a 
second order phase transition for ho/ J = 1.11 and a reen- 
trant behavior of second order appears for ho/ J = 1-2. 
The left panel in Fig. 3a depicts the phase diagram of the 
system in (fcsT c / J — ho/ J) plane with selected values of 
a / J where dotted lines correspond to first order transi- 
tions and solid circles represent the tricritical points. As 
seen from the figure, the system exhibits tricritical points 
and reentrant phenomena for weak randomness and as 
a / J increases then reentrant phenomena and tricritical 
behavior disappear. Furthermore, the critical temper- 
ature of the system decreases and ferromagnetic phase 
region gets narrower with increasing trend in a/J. On 
the left panel in Fig. 3b, we investigate the phase dia- 
gram of the system in (fcsT c / J— a/ J) plane with selected 
values of ho/ J- As seen from the figure, for the values 
of h / J = 0, 0.5, 0.7 and 1.0 the system undergoes a sec- 
ond order phase transition between paramagnetic and 
ferromagnetic phases at a critical temperature which de- 
creases with increasing values of ho/ J. For the values 
of h / J = 1.21,1.22,1.24 and 1.26 the system exhibits 
a reentrant behavior of second order and the transition 
lines exhibit a bulge which gets smaller with increasing 
values of ho/ J which means that ferromagnetic phase re- 
gion gets narrower. Besides, the magnetization profiles 
shown on the right panel in Fig. 3b for ho/ J = 1.21 
with a / J = 0.2,0.4 and 0.6 show that reentrant behav- 
ior disappears with increasing trend in a/J. Random 
field distribution shapes corresponding to the magneti- 
zation profiles in Fig. 3a and Fig. 3b are also illustrated 
in Fig. 3c. 

Finally, in Fig. 4 we show the phase diagram in 
(ho/ J — a / J) plane with constant critical temperature. 
In Fig. 4a, solid and dotted lines correspond to contin- 
uous and reentrant phase transition regions (of second 
and first order), respectively. As we can see from Fig. 
4a reentrant behavior can be observed only for specific 
set of ho/ J and a/J parameters. Fig. 4b is a schematic 
representation of the phase transition characteristics of 
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FIG. 3: (a) Phase diagram (left panel) in (fesT c /J — ho/ J) 
plane with selected values of a/J and variation of magneti- 
zation with temperature (right panel) for selected values of 
ho/ J with a / J — 0.1. Solid and dotted lines correspond to 
second and first order phase transitions and solid circles de- 
note the tricritical points, (b) Phase diagram (left panel) in 
(fcsTc/J — a/J) plane with selected values of ho/ J and vari- 
ation of magnetization with temperature (right panel) with 
specific values of a/J for ho/ J = 1.21. Solid and dotted lines 
correspond to second and first order phase transitions and 
solid circles denote the tricritical points, (c) Double Gaus- 
sian distribution corresponding to the magnetization curves 
depicted on right panels in (b) and (c). 
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the system with double gaussian random field distribu- 
tion. In this figure, the region which is colored in red rep- 
resents ho/ J and a/ J values that cause a reentrant phase 
transition. If we select random field parameters ho/ J and 
a/ J from green colored region the system always under- 
goes only a secon order phase transition and similarly in 
blue region the system always remains in a paramagnetic 
phase for the whole temperature range. Furthermore, al- 
though increasing a/ J value reduces the critical temper- 
ature of the system in general (see Fig. 4a) , it is obvious 
in Fig. 4b that the green region where only second order 
phase transitions take place gets wider for the values of 
(7 < 0.7 and then it gets narrower for a/ J > 0.7. 
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IV. CONCLUSIONS 

In conclusion, we have studied the phase diagrams of 
a spin-1/2 Ising model in a random field by using three 
types of random field distributions. We have adopted 
an effective field approximation that takes into account 
the correlations between different spins in the cluster of 
considered lattice. As a consequence, when the random 
field distribution is single Gaussian the system always 
undergoes a second order phase transition between para- 
magnetic and ferromagnetic phases. This result agrees 
with those of the other works in the literature. For bi- 
modal field distribution, the system exhibits reentrant 
phenomena. On the other hand, our numerical analysis 
clearly indicate that such a field distribution leads to a 
tricritical behavior and this observation contradicts with 
the results of previously published works based on EFT. 
In addition, in the absence of any randomness, numerical 
value of the critical temperature of the system obtained 
within the formalism of the present work agrees well with 
the exact result. Furthermore, double Gaussian form of 
the field distribution as well as single Gaussian counter- 
part has not been examined within the EFT approxi- 
mation before. Hence we have also examined the system 
with this distribution function and we have obtained rich 
phase diagrams displaying reentrant phenomena and tri- 
critical behavior. As a result, we can conclude that all 
of the points mentioned above show that our method 
is superior to the conventional EFT methods based on 
decoupling approximation and therefore the calculation 
results are more accurate. 

We hope that the results obtained in this work may be 
beneficial from both theoretical and experimental points 
of view. 



(a) 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 

o/J 

(b) 

FIG. 4: (a) Phase diagram in (ho/ J — f/J) plane with con- 
stant critical temperature, (b) a schematic representation of 
the phase transition characteristics of the system with double 
gaussian random field distribution. The regions with different 
colors indicates different phase transition characteristics. 
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